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ABSTRACT 

o : 

04 ■ Context. Surprisingly high amounts of H2O have recently been reported in the circumstellar envelope around the M-type asymptotic 

' giant branch star W Hya. This has lead to the speculation that evaporation of icy cometary or planetary bodies might be an elfective 

^ 1' ongoing mechanism in such systems. However, substantial uncertainties remain, as the required radiative transfer modelling is difficult 

, due to high optical depths, sub-thermal excitation and the sensitivity to the combined radiation field from the central star and dust 

. grains. 

' Aims. By performing a detailed radiative transfer analysis, we determine fractional abundances of circumstellar HtO in the envelopes 

around six M-type asymptotic giant branch stars. The models are also used to predict H2O spectral line emission for the upcoming 
Herschel/HIFl mission. 

Methods. We use Infrared Space Observatory Long Wavelength Spectrometer spectra to constrain the circumstellar fractional abun- 
dance distribution of ortho-H^O, using a non-local thermal equilibrium, and non-local, radiative transfer code based on the acceler- 
ated lambda iteration formalism. The mass-loss rates and kinetic temperature structures for the sample stars are determined through 
^ , radiative transfer modelling of CO line emission based on the Monte-Carlo method. The density and temperature profiles of the cir- 

5_( . cumstellar dust grains are determined through spectral energy distribution modelling using the publicly available code Dusty. 

' Results. The determined ortho-HiO abundances lie between 2 x 10""* and 1.5 x 10"^^ relative to Ht, with the exception of WX Psc, 

' which has a much lower estimated ortho-H20 abundance of only 2 x 10"*, possibly indicating HtO adsorption onto dust grains or 

recent mass-loss-rate modulations. The estimated abundances are uncertain by, at best, a factor of a few. 

Conclusions. The high water abundance found for the majority of the sources suggests that either the 'normal' chemical processes 
CO ' are very effective in producing H2O, or else non-local thermal equilibrium atmospheric chemistry, grain surface reactions, or a release 

^ ' of H2O (e.g. from icy bodies like Kuiper belt objects) play a role. However, more detailed information on the physical structure and 

the velocity field of the region where the water vapour lines are formed is required to improve abundance estimates. We provide pre- 
■ dictions for ortho-H20 lines in the spectral window of Herschel/HIFI. These spectrally resolved lines cover a wide range of excitation 

0^ ' conditions and will provide valuable additional information on the physical and chemical properties of the inner stellar wind where 

' H2O is abundant. 
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1 . Introduction for a short while, creates a planetary nebula (Habing[l996|l. The 

lifetime on the AGB is estimated to be a few ~ 10^ yr (Marigo 
Asymptotic giant branch stars (AGB stars) are chemically ^ Girardi i2007i) 
characterised according to their relative abundances of carbon 
^ . and oxygen. Carbon stars have a photospheric C/O ratio > 1, 
■ - ' where as M-t ype ('oxyge n rich ') AGB stars have C/O < 1 (e.g., To date, about 70 circumstellar molecular species have been 
Russel [T934l Beck et al. 11992b . It is beUeved that M-type AGB detected, most of them through transitions at radio wavelengths 
stars evolve to carbon stars by the dredging up of nucleosynthe- (Olofsson l2007l l. H2O is of special importance in the case of M- 
sise d carbo n during the thermally pulsing AGB (TP- AGB) (e.g., type AGB stars, as it is possibly the dominant coolant in the inner 
Iben ll975l Boothroyd & Sackmann ll988. ). The evolution along part of the stellar wind (where it is abundant) with a large number 
the AGB is dominated by the mass loss, starting with a low mass- of far-infrared (FIR) transitions (Truong-Bach et al. 1 19991 1. H2O 
loss rate ~ 10"** Mq yr"' and ending in an intense wind with rates line emission is also a potentially important probe of the physi- 
up to lO"'* Mq yr"' . The photosphere and expanding circumstel- cal and chemical properties of the inner regions of CSEs. H2O is 
lar envelope (CSE) are effective producers of a variety of molec- the main reservoir of circumstellar O, and whereas CO is abun- 
ular species and microscopic d ust par ticles, providing up to 80% dant in all AGB stars, chemical models predict H2O to be one of 
of the dust in galaxies (Whittet[T952|l. When the star reaches the the most abundant species only in O-rich stars (e.g., Cherchneff 
end of the AGB, most of the initial mass has been lost, leaving 2006). However, H2O has also been detected in the nearby C- 
behind only a C/O core that ionises the suiTounding material and, ric h AG B star IRC+10216 (Melnick et al. 120011 Hasegawa et 
al. l2006l l. but here the estimated abundance is comparatively low, 
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The atmosphere of the Earth is mostly opaque at H2O Hne 
wavelengths and it is therefore difficult to observe H2O from the 
ground, in particular in its ground vibrational state. However, 
Menten et al. ( 12006b managed to observe two vibrationally ex- 
cited lines in the V2 = 1 bending mode (523 - 616 at 336.2 GHz 
and 651 - 752 at 293.6 GHz) towards the red supergiant VY 
CMa. H2O masers from the ground and excited vibrational states 
can also be detected from the ground (Menten & Melnick ,1989D . 
Observations from the Infrared Space Observatory (ISO) are 
not hindered by the atmosphere, and the Short- and Long- 
Wavelength Spectrometer (S WS and LWS) data show a rich H2O 
spectrum in M-type AGB stars (e.g., Truong-Bach et al. !1999! 
Barlow et al. .1996 , Neufeld et al. 119961 Justtanont et al. |2005D . 
The ortho-H20 ground state hne at 557 GHz was also observed 
by the SW AS (Melnick et al. l2000l l and Odin satelUtes (Nordh et 
al. l2003l l. In contrast to the CO radio lines, the H2O line emis- 
sion originates from closer to the star, providing information on 
the warm, high-density innermost layers of the CSE (Truong- 
Bach et al. 11999b . However, the high optical depths and the sub- 
thermal excitation of H2O in CSEs make it a challenge to accu- 
rately model the observed H2O spectral lines (e.g., Justtanont et 
al. l2505l l. 

We present here the results of radiative transfer modelling 
of the ortho-H20 line emission in six M-type AGB stars with 
mass-loss rates between 1 x lO^^Moyr"' and4x lO^^'Moyr"'. 
Detailed radiative transfer models are created using the acceler- 
ated lambda iteration (ALI) method (Rybicki & Hummer 1 1 99 1 1 
11992b . The model results are fit to LWS ISO observations be- 
tween 43 and 197 //m. The mass-loss rate and kinetic tem- 
perature structure are determined through CO radio line mod- 
elling, and the dust temperature and dust density structure are 
determined using Dusty (Ivezic et al. 11999b . The results show 
the ability of ALI to model circumstellar H2O lines. We also 
make predictions for H2O lines in the frequency range of the 
upcoming Herschel/HIFI mission. These observations will in- 
clude information on the line profiles and can be used to set 
firmer constraints on the physical structure of the inner CSE 
and the H2O abundance. Previous attempts to determine the wa- 
ter vapour content in stars included in our sample have been 
made for W Hy a (Bar low et al. [T9961 Zubko & Elitzu r l2000l 
Justtanont et al. l2005l l and R C as (Tru ong-Bach et al. [T999] l. 
Gonzalez- Alfonso & Cernicharo ( 11999b used the large-velocity- 
gradient (LVG) method to model the 183 GHz maser line in 
O-rich evolved stars, to study the line intensity dependence on 
source properties and physical conditions. Ryde & Eriksson 
( 12002b model the 2.6 - 3.6 jjm spectrum for R Dor observed with 
the S WS on ISO, showing the dominance of water vapour in the 
spectrum. 

Determining the para-H20 abundance would be a straight- 
forward task. However, considering the uncertainty in the abun- 
dance estimates, this would not give any firm constraints on the 
ortho-to-para ratio or the total H2O abundance. 

In Sect. 12] the ISO observations are briefly described. 
Section [3] contains a description of the dust and CO emission 
line modelling. In Sect. |4]we describe the H2O line model, and 
in Sects. |5]and|6]we present and discuss the results, respectively. 
Appendix |A] contains a description of the accelerated lambda 
technique, adopted in our radiative transfer model, in more detail 
and benchmark tests are presented. All abundances mentioned 
refer to the fractional abundance compared to H2, unless stated 
otherwise. 



Table 1. ISO observations for the 6 M-type AGB stars. 



Source 



AOT 



Observer ID 



Date 



Int time (sec) 



TXCam 

RCas 

RDor 

WHya 

WX Psc 

IKTau 



LWSOl 
LWSOl 
LWSOl 
LWSOl 
LWSOl 
LWSOl 



MBARLOW 
MBARLOW 
KERIKSSO 
MBARLOW 
MBARLOW 
MBARLOW 



1997-10-10 
1997-06-06 
1997-07-01 
1997-08-02 
1997-06-15 
1997-09-02 



3428 
2204 
1330 
1778 
2796 
3430 



2. Observations 

We have used archived ISO LWS spectra (43 - 197 jum) for six 
M-type AGB stars (Table[T]) to identify and measure the intensity 
of ortho-H20 spectral lines in the ground vibrational state. For 
all observations the LWSOl AOT was used (Clegg et al. 11996b . 
The spectra were sampled at 1/4 of a resolution element, one el- 
ement being 0.3 jjm in second order {A < 93 /im) and 0.6 ;um in 
first order (A > SOyum). The analysis of the data was done us- 
ing the ISO Spectral Analysis Package (ISAP). The data were 
cleaned from cosmic rays by excluding deviant points that were 
only present in individual scans. All scans for each detector 
were median averaged. In order to get a straight (continuum- 
subtracted) spectrum, a baseline (typically of order 2) was fit to 
the data and subtracted. Finally, the ten detectors were combined 
to form a complete spectrum. A section of the merged spectrum 
for R Dor with labelled ortho- and para-H20 transitions is shown 
in Fig. [U Table |2] gives the measured ortho-H20 line fluxes. 
The observed fluxes are determined by fitting a Gaussian pro- 
file to the observed ortho-H20 lines using available tools within 
ISAP. The error in the baseline fit dominates over the error in 
the Gaussian fit, and is ~ 30% on average (in terms of the line 
intensity). The fit is more difficult in the short-wavelength part of 
the spectrum due to the high noise levels. The measured values 
lie within + 20% of those of previous studies of H2O emission 
using ISO spectra in W Hya (Barlow et al. 1 1 996b and R Cas 
(Truong-Bach et al. [T999b . 



3. Circumstellar dust and CO line modelling 

3.1. Dust emission modelling using Dusty 

The dust optical depth and dust temperature profile were deter- 
mined using the continuum radiative transfer code Dusty (Ivezic 
et al. 119991 ). The observational constraints are in the form of 
spectral energy distributions (SEDs) and consist of 2MASS and 
IRAS fluxes. A 'standard model' containing amorphous silicate 
dust was used, with optical constants adopted from Justtanont 
& Tielens (119921) . The standard model assumes the mass loss to 
be spherically symmetric with a constant rate and expansion ve- 
locity, resulting in a density profile with p oc r^^. Prompt dust 
formation is assumed at the condensation temperature Tc, which 
is treated as a free parameter. Additionally, the dust grains are as- 
sumed to be spherical with a radius of 0. 1 fim. The dust optical 
depth and dust condensation temperature were adjusted in order 
to find the best fit to the SED using a ;^'^-statistic. A more de- 
tailed description of the modelling procedure is given in Schoier 
et al. ( [2002b . 

The luminosities of the Miras are derived using a period- 
luminosity (PL) relationship based on bolometric magnitudes 
(Feast et al. 11989b . The distances are obtained from the appar- 
ent bolometric magnitude obtained in the SED modelling. For 
semi-regular objects, Knapp et al. ( 12003b give a PL relationship 
based on K-band magnitudes. Since this only uses information in 
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Fig. 1. The ISO LWS spectram of R Dor between 95 and 185 jum, including the position of ortho-H20, para-H20 lines, and CO 
lines (transitions J - 15 - 14 to 7 = 27 - 26) . 

Table 2. Measured and modelled ISO line fluxes of ortho-H20 [in units of 10 ""W cm Lines without fluxes were either too weak 
to detect or otherwise unusable, due to blending with other lines or unacceptable line fits. The dominant error in the observed fluxes 
is due to the baseline subtraction and is ~ 30%. 
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one point of the SED (the K-band magnitude), we instead used 
the revised Hipparcos distances for W Hya and R Dor (Knapp et 
al. |2003]), and determined the luminosity through the SED mod- 
elling. 



3.2. Modelling of CO lines 

The CO line modelling is based on the Monte-Carlo method pre- 
sented in detail in Schoier & Olofsson ( 12001b and Olofsson et al. 
(120021) . The radiative transfer model includes the lowest 40 ro- 
tational levels in the ground and first vibrational states. The col- 
lision rates are taken from Flower ( 120011 ) for CO-H2, and they 
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are extrapolated to higher rotational temperatures as was done 
in Schoier et al. (I2005I I. The collisional rates between CO and 
ortho- and para-H2 were weighted assuming an ortho- to para- 
H2 ratio of 3. A CO fractional abundance (CO/H2) of 2 x 10 is 
adopted, a value typically used for O-rich AGB stars. The model 
includes cooling terms from the adiabatic expansion of the en- 
velope and radiative line cooling through CO and H2. The dom- 
inant heating term is given by dust-gas collisions. The density 
distribution and temperature structure of the dust grains are taken 
from the dust radiative transfer modelling (see Sect. 13.1b . while 
the parameters for grain size, grain mass density and dust-to-gas 
ratio are combined into one parameter h. The /i-parameter is, in 
principal, a free parameter in the CO modelling, and is chosen to 
be consistent with previous models (Olofsson et al. 120021 W Hya 
and R Dor; Gonzalez Delgado et al. 120031 R Cas; Ramstedt et 
al. 120071 TX Cam, WX Psc, and IK Tau). The o uter ra dius is 
set by the photodissociation of CO (Stanek et al. '1995 , based 
on results by Mamon et al. |1988t . The expansion velocity and 
micro-turbulent velocity (0.5 km s ' doppler width), with local 
thermal (Gaussian) line broadening, are also set in the model. 
This leaves the mass-loss rate as the remaining free parameter 
that affects the model line fluxes. The best fit model is then de- 
termined by fitting the model line fluxes to the observed fluxes 
and minimising where 

7 \ 1 (^obs.i ~ ^mod.i) 

err 
1 

and /obs and I^od are the observed and modelled line fluxes, re- 
spectively, CT; the uncertainty in observation /, and the sum goes 
over all lines included. 

The line profile is not taken into consideration, but we note 
that in the case of W Hya, Justtanont et al. ( I2005I I had to reduce 
the CO photodissociation radius and increase the mass-loss rate 
in order to fit the line profile. However, there is sufficient uncer- 
tainty in the underlying assumptions of the circumstellar model 
that we prefer to treat all objects in the same way and avoid in- 
dividual adjustments of, e.g., the CO envelope size. 

The uncertainty in the mass-loss rate is estimated to be 
+ 50% within the adopted model (Schoier & Olofsson 12001b . 
Previous mass-loss-rate estimates by Gonzalez Delgado et 
al. (12003b and Olofsson et al. (12002b differ somewhat from the 
results in Table |4] The difference is less than 10% for the low- 
mass-loss-rate and intermediate-mass-loss-rate objects, but more 
for the high-mass-loss-rate objects WX Psc and IK Tau (+ 50%). 
This is due to the combination of new adopted distances, inclu- 
sion of new data, updated collisional rate coefficients and the 
inclusion of thermal radiation from dust grains in our CO mod- 
els. The observed CO line fluxes, used in the model fitting, are 
presented in Table |3] IK Tau, WX Psc and TX Cam were mod- 
elled by Ramstedt et al. (12007b . The parameters derived from the 
dust and CO modelling are presented in Table |4] 

4. Circumstellar H2O line modelling (ALI) 

We have developed a detailed radiative transfer code based on 
the accelerated lambda iteration method, in order to accurately 
model circumstellar H2O line emission (Bergman, internal re- 
port). The Monte Carlo method is not well suited for such prob- 
lems due to its slow convergence at high optical depths. The ALI 
code was previously used in Justtanont et al. (12005b . but a num- 
ber of modifications have subsequently been made to be able to 
treat the inclusion of dust emission and vibrationally excited en- 
ergy levels. The accelerated lambda technique and its numerical 



Table 3. Observed CO line fluxes in main beam brightness scale. 
Observations were taken at Onsala Space Observatory, JCMT, 
and SEST (see references for details). 



Source 


y=i->o 

[Kkms-'l 


y = 2^ 1 

[Kkms-'] 


7 = 3^2 
[Kkms-'] 


7 = 4^3 
[Kkms-'l 


TXCam 


20.1'' 


60.6^ 


71. r 


148.7'' 


RCas 


8.4" 


32.1' 


100.2' 


89.6* 


RDor 


5.0" 


41.6" 


70.4" 




WHya 


0.83" 


18.7'' 


40.2'' 


44.6'' 


WXPsc 


44.4^ 


80.6" 


70. r 


83.8" 


IK Tau 


46.2^ 


95.r 


124.5'' 


129.6" 



Kerschbaum & Olofsson i ri999ll , Delgado et al. ( |2003]| , " JCMT 
archive, Justtanont et al. l i2005t . = Ramstedt et al. ( I2007t 



implementation is described in more detail in appendix lAl in- 
cluding results of benchmarking tests against the standard prob- 
lems presented in van Zadelhoff et al. ( 12002b . 

4. 1 . Dependence on numerical parameters 

ALI calculates the mean intensity at each radial point in the CSE, 
including the contribution from all other points in the envelope. 
This is done by paying special attention to a proper sampling 
of the central source. Typically, 50 radial points (distributed in 
a logarithmic way) and 32 angular points were used. The de- 
pendence on the radial and angular sampling was tested for the 
low mass-loss-rate object R Dor and intermediate mass-loss-rate 
object TX Cam by increasing the two parameters by a factor of 
two. The differences in the modelled line fluxes were small in 
both cases, less than 2%. 

After each iteration the new level populations are calculated, 
using a user-defined fraction of the old level populations (here 
set to 0.8) to prevent the models from diverging too easily. The 
models are said to converge when the average change in the level 
populations between iterations is less than a defined limit. This 
limit was also tested by running models with different limits for 
R Dor and IK Tau. For an average change between iterations in 
the level populations of less than 10"^ the models changed by 
less than 0.1% This limit was subsequently used for all models. 

4.2. Modelling of RjO lines 

The radiative transfer model includes the 45 lowest levels in the 
ground and first excited vibrational states (the bending mode 
V2 - I, 6.3 /im or «2300K above the ground state) of ortho- 
H2O, while excitation through the stretching modes can be ne- 
glected. For the high mass-loss-rate objects IK Tau and WX Psc, 
only the ground vibrational state was included in order to im- 
prove the convergence. The effect of excluding the V2 state on the 
resulting model line fluxes was tested and found to be insignifi- 
cant (< 4%) for these high-mass-loss-rate objects (see Sect. 15.3b . 
The collisional rates in the ground state are taken from the rates 
of H2O with He, corrected by a factor of 1 .4 to account for col- 
lisions with H2 (Green et al. 1993). Phillips et al. ( 11996b suggest 
that this may lead to rather different rates compared to those ob- 
tained when including H2 in the calculation. However, they only 
calculate collisional rates between H2O and H2 up to a temper- 
ature of 140 K, while typical temperatures in the CSEs of AGB 
stars can be up to 1000 K. The importance of collisions within 
the V2 state and between the ground and V2 state was also tested 
(see Sect. 15.31 ). Rotational collision rates within the first excited 
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Fig. 2. Modelled ISO ortho-H20 lines for TX Cam. The dot- 
dashed line marks the stellar vlsr, the solid line shows vlsr ± 
Vexp- Self-absorption and weak P-Cygni profiles are present in 
the model line profiles. 



state are taken to be the same as for the ground state, while col- 
lisions between the ground and excited V2 state are based on 
the ground state rotational collision rate coefficients scaled by 
a factor of 0.01 . The molecular data for H2O was taken from the 
Leiden Atomic and Molecular Database^ (LAMDA) (Schoier et 
al. 120051) . 

As with the CO models, the density distribution and temper- 
ature structure of the dust grains are taken from the dust radiative 
transfer modelling (see Sect. 13. 11 1. From the CO line modelling 
the mass-loss rate, expansion velocity, and kinetic temperature 
profile are included (see Sect. l3.2l and TablelDi. The inner radius 
of the CSE is set to the inner radius of the dust envelope, i.e. the 
dust condensation radius. Setting a smaller radius does not affect 
the line fluxes. For all models, a Gaussian distribution of circum- 
stellar H2O is assumed. The maximum radius is defined by the 
e-folding radius (i.e., where the abundance has decreased to 
37%) and is set to 3 x rg. By varying r^, the size of the circum- 
stellar H2O envelope is treated as a free parameter 

An estimate of re is given by the results of theoretical models 
of H2O photodissociation (Netzer &l Knapp [T987| l, where the ra- 
dius at which the abundance of OH (the photodissociation prod- 
uct of H2O) is highest is given by 

The remaining parameter to fit is the ortho-H20 abundance 
/0-H2O relative to H2. The model lines are fit to the line fluxes 
of the observed spectra in the same way as for the CO emission 
line models. The best fit model was determined, as for the CO 
modelling, by minimising as defined in Eq. [T] with and 
/0-H2O as free parameters. 



' available at |http://www.strw.leidenuniv.nl/~moIdata| 



5. Results 

5.1. Modelled lines and ortho-H20 abundances 

For each object a grid, of models with varying ortho-H20 abun- 
dance and H2O envelope size was calculated. The resulting 
maps are shown in Fig.|3] Between 4 and 16 lines are matched 
for each object. The fluxes of the best-fit modelled ortho-H20 
lines are shown in Table |2] together with the observed fluxes. 
For the low-mass-loss-rate objects the abundance and rg can be 
reasonably well determined, while setting constraints for r^. be- 
comes increasingly diflicult with increasing mass-loss-rate. In 
the high-mass-loss-rate objects, the emitting region is excitation 
limited. Consequently, the emission is not very sensitive to the 
envelope size. As a result, it is difficult to constrain the size of 
the of the H2O envelope for these objects. For the low-mass-loss- 
rate objects, the excitation can be high throughout the envelope, 
thus making the emission sensitive to the envelope size, i.e. the 
emission can be said to be photodissociation limited. The grids 
are terminated in envelope size when an increase in does not 
lead to a significant change in line fluxes anymore, or when the 
lo" contour closes. For IK Tau and WX Psc a lower limit could 
not be determined, as smaller envelope sizes require a higher 
abundance, for which the models could not be made to converge. 
The;^^^ maps are truncated at 0.1 x tnk for these objects. In 
the cases where an outer radius could be determined, the radius 
given by Eq.|2]lies within the 1 cr contour of the lowest x^ model 
radius. Table|4]gives the radius from the theoretical model (rNK), 
the radius of the lowest x^ model (rg), and the corresponding 
abundance estimates (/nk and /g, respectively). From here on, 
the 'best-fit' models refer to the results obtained using tnk from 
Eq.E] 

The resulting abundances of circumstellar ortho-H20, from 
fitting the model lines to the observed ISO lines, lie between 
2 X 10"* and 1.5 x 10"^^ relative to H2 (Table Si. Four of the 
modelled ISO lines of TX Cam are shown in Fig. |2]as an ex- 
ample of an intermediate-mass-loss-rate object. Self absorption 
is present in all the lines, whereas weak P-Cygni profiles are 
seen in the higher-frequency lines, where the continuum emis- 
sion from thermal dust grains is stronger. Since the line shapes 
are not spectrally resolved in the ISO spectra, only the fluxes 
were fit. HIFI will be able to resolve the line profiles, and these 
can then be taken into consideration when fitting the models 
(Sect.l531l. 

Figure |4] shows the kinetic temperature profile compared to 
the excitation temperature for four lines for W Hya and IK Tau 
(upper panels), and the tangential optical depth at line centre for 
the same lines as a function of distance from the central star 
(lower panels). It should be pointed out that the optical depth in 
the line wings usually is considerably higher Two of the lines 
are observed by ISO, while three of the lines will be observ- 
able by HIFI. In general, the Unes are subthermally excited, and 
the Jk_k^ - Jk-K^ - 532 - 441 transition at 620.9 GHz shows 
maser action. The lines are more thermalised in the high-mass- 
loss-rate object IK Tau. The tangential optical depth shows the 
regions where the lines are created, and this indicates that the 
ISO observations do indeed probe the inner regions of the H2O 
envelope. W Hya has a smaller envelope due to the lower mass- 
loss-rate and consequently the lines are formed closer to the star. 

5.2. Constraining the H2O abundance distribution 

Several test models were calculated in order to see how the re- 
sults depend on the ortho-H20 abundance, mass-loss-rate, and 
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¥ig.3.x^ maps in ortho-H20 abundance and radius of the H2O envelope for all stars in the sample. The contours give the 1, 2 and 
3cr limits. The stars mark the best model based on the theoretical radius given by Eq.|2] the dots mark the models with the lowest 
values. The;^'^^^ values in the figures refer to the best model fits. 



Table 4. Results of the dust emission, and CO and H2O line modelling (see text for details), tnk and r^. are the theoretical model 
radius and the e-folding radius, respectively. The corresponding H2O abundances for both radii are given under /nK and /e, 
respectively. The;^'^^^^ values refer to the best H2O model fits using tnk. Likewise, diff gives the average difference of the modelled 
H2O line fluxes using tnk to the observations in %. 



Source 


Type 


L" 






M 


^exp 


'■nk 




/nk 


/e 




diff 






[Lo] 


[pc] 


10/jm 


[10-«MoW-'] 


[kms-'] 


[lO'^cm] 


[lO'^cm] 


[10-"] 


[10--*] 




% 


TXCam 


M 


11900 


440 


0.4 


7.0 


18.5 


13.1 


20.0 


3.0 


3.0 


0.6 


-6.8 


RCas 


M 


8725 


172 


0.05 


0.9 


10.5 


3.9 


3.0 


3.0 


3.5 


0.3 


-4.8 


RDor 


SRb 


6500 


59 


0.03 


0.2 


6.0 


1.7 


1.1 


2.0 


3.0 


0.7 


-1.9 


WHya 


SRa 


5400 


78 


0.07 


0.1 


7.2 


1.1 


2.2 


15.0 


15.0 


0.6 


-10.4 


WX Psc 


M 


14600 


720 


3.0 


40 


19.3 


43.0 




0.02 




2.7 


-7.3 


IKTau 


M 


10400 


300 


1.2 


10 


19.0 


17.0 




3.5 




1.2 


-8.5 



" From period-luminosity relation (Feast et al. ll989t , except for R Dor and W Hya, where the Hipparcos distances and the dust modelling are 
used. 

Determined in the dust modelling using L, except for R Dor and W Hya, where the Hipparcos distances are used (Knapp et al. l2003t . 



size of the circumstellar envelope. This was done by varying one 
of the parameters, while all other parameters were held constant. 
Table |5] shows the results in the cases of R Dor and TX Cam, as 
examples of low and intermediate-mass-loss-rate objects. 

Increasing and decreasing the abundance by a factor of two 
changes the line fluxes on average by +20% and ±10% for the 
low- and high-mass-loss-rate objects, respectively (i.e. less than 
the observational error in the ISO lines). In general, the error in 
the H2O abundances is estimated to be of the order +50% within 
the adopted model. 

Varying the mass-loss rate by +50% changes the line fluxes 
by a; + 30%. The change in mass-loss rate aff'ects the density dis- 
tribution, and hence the efficiency of the collisional excitation, 
making the lines more sensitive for changes in the mass-loss rate 
than for changes in the abundance. 

Increasing and decreasing rg by a factor of two changes the 
line fluxes by less than ±20%. The emitting region is excitation 
limited, instead of limited by the size of the H2O envelope, mak- 
ing it difficult to provide good constraints on re . The position of 



OH masers for W Hya (0.8 - 1 .4 x lO'^ cm from the central star, 
see Szymczak et al. (fTg98] l) and for R Cas (1.3 - 2.6 x lO'^ cm 
from the central star, see Chapman et al. ( I1994I I). indicate that 
the H2O envelope might actually be smaller than given by Eq.|2] 
We conclude that the saturation of the lines, and hence the in- 
sensitivity of the integrated intensity to the abundance, is the ma- 
jor source of uncertainty in the determined abundances. Taking 
into account the uncertainty in the mass-loss rate (Ramstedt et 
al. 2007), the line fluxes, and in the model-dependent assump- 
tions (such as the molecular distribution and velocity field), the 
determined abundances are accurate to within a factor of approx- 
imately 5. 

5.3. Dependences on 'fixed' parameters 

Additional test models were calculated to examine the depen- 
dence of the line fluxes on the 'fixed' parameters (see Table|5]l. 
As before, only one parameter was varied, while all other pa- 
rameters were held fixed. The effect of excitation to the first 
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10^* 2x10^* 5x10^* 10^^ 
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Fig. 4. Kinetic temperature profile and excitation temperatures for four modelled lines (indicated by their Jk_k+ - Jk-K^ quantum 
numbers) in W Hya and IK Tau (upper panels), and the tangential optical depths (at line centre) for the same lines (lower panels). 
The lines include a high transition line observed with ISO (523 - 4i4 at 45 jum), a low lying transition observed with ISO (3o3 - 2 12 
at 174.6 yum), and the 557 GHz line (lio - loi) observed by Odin and SWAS (the latter two can also be observed by HIFI). A 
line, exhibiting maser action, which can be observed by HIFI is shown (532 - 441 at 621 GHz); the dashed line indicates the region 
where the level populations are inverted, and shows the absolute value of the optical depth in the lower panel. For radii larger than 
3 X lO'^cm, numerical noise dominates the 532 - 441 transition in IK Tau and the model is consequently terminated at this radius. 



excited vibrational state was tested for the four stars for which 
models including this state converged. Excluding excitation to 
this state leads to a stronger decrease of the line fluxes for the 
low-mass-loss rate objects than for the intermediate and high- 
mass-loss rate objects. Hence the error introduced by not includ- 
ing the first excited vibrational state for IK Tau and WX Psc is 
expected to be well below the uncertainty due to the observed 
line fluxes. For the low-mass-loss-rate objects, the inclusion of 
vibrational excitation results in significantly lower abundances 
compared to models that do not include the excited vibrational 
state (Justtanont et al. 12005b . The variation of the line fluxes 
when excluding the dust radiation field and the central star was 
also examined in test models for the low-mass-loss-rate object 
R Dor and intermediate-mass-loss-rate object TX Cam (Table|5]). 
It must be kept in mind that the changes in line fluxes all lie 
within (or close to) the observational uncertainties. This is es- 
pecially true for TX Cam, where the change between different 
models is at most 12% on average. 

Including the rotational collision rates within the V2 = 1 state 
and between the ground state and first excited state, as described 
in Sect. 14.21 does not change the fluxes significantly for any 
of the stars, showing that the dust radiation field and radiation 
from the central star are more important for the excitation of the 
V2 = 1 state than collisions. The role of collisions within the 



vibrational ground state is tested by setting all collisional coef- 
ficients to zero, decreasing the observed line fluxes, on average, 
by 17% in the case of R Dor and 58% for TX Cam. As the lines 
are generally subthermally excited, the error in the collision rates 
in the ground state is likely to affect the modelled line fluxes (see 
Sect.feb. 

Since the ISO H2O emission lines probe the inner parts of 
the CSE, the expansion velocity in the region where the lines are 
created may not have reached its terminal value as determined 
in the CO line modelling (see Sect. l3.2b . [Note that changing the 
expansion velocity, while keeping the mass-loss-rate constant, 
effectively changes the density profile.] Decreasing the expan- 
sion velocity by 20% only changes the line fluxes by -h6% for 
both TX Cam and R Dor, showing that the uncertainty in the ex- 
pansion velocity does not contribute significantly to the error in 
the estimated abundances. 



5.4. H2O line cooling 



In addition to CO and H2 cooling by line emission, H2O could 
also be an important coolant in the inner part of CSEs where it is 
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Table 5. Parameter dependence of the model ISO line fluxes for 
the low-mass-loss-rate object R Dor and the intermediate-mass- 
loss-rate object TX Cam. Column 2 gives the change of the vari- 
ous parameters relative to their best-fit value (Table|4]i. Columns 
3 and 4 give the average difference of all model line fluxes com- 
pared to the best-fit model. 



param. 


change 


RDor 


TX Cam 


,/e 


+ 100% 


+20% 


+ 10% 




-50% 


-20% 


-10% 


^mod 


+ 100% 


+20% 


+20% 




-50% 


-20% 


-20% 


M 


+50% 


+30% 


+30% 




-50% 


-30% 


-30% 




-20% 


+6% 


+6% 


T'kin 


±20% 


±15% 


±15% 


collisions 


excl. 


-17% 


-58% 


V2 


excl. 


-35% 


-A% 


dust 


excl. 


-11% 


-12% 


star 


excl. 


-19% 


-3% 



abundant. From the excitation analysis, the cooling (or heating) 
arising from line emission can be estimated using 

A = wh, ^ ^(nm« - nuYuiViVui, (3) 

where «„ and «/ are the level populations in the upper and lower 
levels participating in the transition (with rest frequency v,,/) and 
■y„/ and y;„ are the collisional rate coefficients (van der Tak et 
al. l2007l l. The cooling rate A in ergs ' cm is defined as pos- 
itive for net cooling. In Fig. |5]the line cooling from 0-H2O is 
compared with that from CO and H2 for the W Hya and IK Tau 
models used in the present analysis. In the region where H2O 
is abundant its cooling rate will be more than an order of mag- 
nitude higher than that from CO. Due to the constraints set by 
the CO-line intensity ratios, the kinetic temperature profile de- 
termined through the CO line modelling is most accurate in the 
region in which the CO lines are emitted. Any major cooling 
due to H2O must therefore be partly compensated by additional 
heating, by adjusting the /i-parameter so that the kinetic temper- 
ature structure is consistent with the observed CO rotational line 
intensities. 

To include molecular cooling from water self-consistently, 
the energy balance equations for CO and H2O must be solved si- 
multaneously. For a complete description of the energy balance, 
other processes (such as a complicated heating term) would have 
to be included as well, which is beyond the scope of this paper. 
A crude estimate of how H2O line cooling and other processes 
might affect the modelled line intensities is done by changing the 
temperature profile resulting from the CO modelling by + 20%. 
This changes the line intensities on average by less than +15%, 
without changing the line ratios. 

5.5. Predictions for HIFI 

HIFI (Heterodyne instrument for the far infrared) is a hetero- 
dyne receiver onboard the Herschel space telescope, scheduled 
for launch in 2008. The instrument is a high-resolution (up to 
10^, 0.03 - 300 km s ') heterodyne instrument designed to pro- 
vide a wide and continuous frequency coverage. The frequency 
range covered includes a large amount of molecular lines (de 
Graauw et al. 120051 ). We have modelled 19 ortho-HaO lines cov- 
ering a range of excitation and optical depths within the range 
557 GHz to 1885 GHz. 




r [cm] r [cm] 



Fig. 5. The emissivity A in ergs"' cm"^, scaled by (r/10''*cm)"', 
of the line cooling rate from 0-H2O, CO and H2 in the circum- 
stellar envelopes around W Hya and IK Tau. There is a net cool- 
ing (A > 0) throughout the envelopes for these molecules, except 
very close to the star in the case of CO in W Hya (A < indicated 
by a dashed line). 



The sensitivity limit of HIFI for a 5 cr line detection after 1 
hour of integration at a resolution of R = 100000 is between 2.2 
and22.0x 10-^°Wcm-2Hz-i (de Graauw et al.|2005ll, depend- 
ing on the frequency band. The modelled HIFI lines generally 
have fluxes of 10"^^ W cm ^ Hz 'or higher, making them ob- 
servable without difficulty. Typically, between 1 and 10 hours of 
integration time is required to obtain a sufficient signal-to-noise 
ratio, even for weaker lines in the first excited vibrational state. 

Figure |7] shows a selection of the modelled lines obtained 
from the best-fit models for TX Cam, W Hya and IK Tau as 
given in Table |4] (based on the theoretical H2O envelope radius 
given in Eq. |2]i. The dot-dashed line marks the stellar vlsr and 
the solid line marks Maser action is present in the 532 -4-4i 

and 634-541 transitions at 621 GHz and 1 158 GHz, respectively, 
in all three sources. 

The line profiles for models with lower mass-loss rates gen- 
erally appear to suffer more from self-absorption and show sig- 
nificant emission outside the range visR+Vexp- The CSEs of these 
sources are relatively small and the optical depths at the outer 
parts of the CSEs are around, or larger than, unity (see Fig. |4]i. 
This leads to effective self-absorption on the blue-shifted side, 
but also to a contribution of the local line profiles to the line 
wings. This is particularly pronounced in the lines where the op- 
tical depths at line centre can reach up to ~100 (Fig.|4li for these 
kinds of objects, such as lio - loi and 3o3 - 2i2 (Fig. |7]i. The 
broadening of the line profiles is further enhanced by the fact 
that the thermal motion, v,/,, is of the same order or larger when 
compared with the microturbulent motion, Vturb, in the emitting 
region. The result is that optically thin lines can also have signifi- 
cant line broadening (Fig.|6]l. These effects are less significant in 
the high-mass-loss-rate objects, as the optical depths in the outer 
parts of these CSEs, which contribute most to the observed emis- 
sion, are around, or below, one. Also, for the high-mass-loss-rate 
objects v,i, « v,i,rb in the outer parts of the envelope. 

The sensitivity of the HIFI lines to changes in various pa- 
rameters has also been tested. A change in abundance by +50% 
changes the HIFI line fluxes on average by x +30% compared 
to the best-fit model. This compares to the 10 - 20% change 
in the ISO lines, indicating that the lines in the range of HIFI 
are somewhat more sensitive to the H2O abundance (see below). 
Increasing and decreasing rg by a factor of two change the line 
fluxes on average by less than +16%. The line fluxes for the HIFI 
lines are more sensitive for a -20% change of the expansion ve- 
locity than the ISO lines (+20% on average; however, optically 
thin lines are more sensitive, see below). Resolved line profiles 
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Table 6. Model maximum intensities for 3 optically thin lines 
within the range of HIFI (in units of 10"^" W cm"^ Hz '). 



Trans 


V 


RDor 


TX Cam 




[GHz] 






854 - 761 


1169 


10 


4 


854 - 927 


1596 


3 


1 


845 - 752 


1885 


320 


40 



Table 7. Parameter dependence of the model fluxes of three op- 
tically thin lines observable by HIFI for the low-mass-loss-rate 
object R Dor. Column 2 gives the change of the various pa- 
rameters relative to their best-fit value (Table|4|. The remaining 
columns give the change of the model fluxes compared to the 
best-fit fluxes for the respective lines. 



param. 


change 


854 - 761 


854 - 927 


845 - 752 


/e 


±50% 


±30% 


±6% 


±55% 


D 

"mod 


+ 100% 


+4% 


+5% 


+2% 




-50% 


-2% 


-2% 


-A% 


M 


+50% 


+77% 


+ 10% 


+ 153% 




-50% 


-25% 


-7% 


-45% 


^exp 


-20% 


+28% 


+5% 


+55% 




+20% 


+21% 


+ 1% 


+55% 




-20% 


+9% 


-1% 


+20% 



will be crucial for getting information on the dynamical structure 
of the inner CSE. 

Of the 19 modelled lines, 5 are optically thin throughout the 
CSE. Two of these are maser lines (the 532-441 and 634-541 tran- 
sitions at 621 GHz and 11 58 GHz, respectively). The remaining 
three lines are shown in Table |6] with their predicted maximum 
intensities. The 5 cr line detection limit after 1 hour of integra- 
tion time at = 100000 is a; 2 x lO"^** W cm^- Hz^' , making it 
difficult to detect some of these lines. 

The sensitivity of these lines to the ortho-HaO abundance, 
the outer radius, mass-loss rate, kinetic temperature profile, and 
expansion velocity was tested for R Dor and TX Cam. The re- 
sulting change in line fluxes is presented in Tables |7] and [8] for 
R Dor and TX Cam, respectively. Except for the 854 - 92? tran- 
sition, the optically thin lines are much more sensitive to the 
abundance, mass-loss rate and expansion velocity than the ISO 
lines are. The P-Cygni profile causes the insensitivity of the inte- 
grated flux in the 854 - 927 transition. The high sensitivity to the 
expansion velocity is due to the change in the density profile as 
a result of keeping the mass-loss rate constant while decreasing 
the expansion velocity. They are only slightly more sensitive to 
a change in the kinetic temperature profile compared to the ISO 
lines, and not at all sensitive to the outer radius. For R Dor, tran- 
sitions 321 -3 12, 2 12 - 1 01, and 221 -2i2 at 1163 GHz, 1670 GHz, 
and 1661 GHz, respectively, are more sensitive to the change of 
the outer radius (with up to a 30% change in line flux), but the 
emission region is essentially excitation limited. 

An example of how the 845 - 752 line profile changes be- 
tween the different models for R Dor (top) and TX Cam (bot- 
tom) is shown in Fig. |6] The left panels show models in which 
the parameters were reduced as indicated in Tables |7] and [8] and 
the right panels show the models with increased parameters. The 
line profiles change mainly in intensity. However, for TX Cam 
the optical depth in the 845 - 752 transition is just below one 
in the best fit model. The line therefore changes from optically 
thick to optically thin for the different models, and as the optical 
depth decreases, the line becomes increasingly flat-topped. The 



Table 8. Same as Table|2l but for the intermediate-mass-loss-rate 
object TX Cam. 



param. change 854 - 76i 854 - 927 845 - 752 

±50% ±30% ±7% TSOW 

i^mod +100% +20% +3% +6% 

-50% -14% -2% -6% 

M +50% +70% +14% +71% 

-50% ^4% -11% -53% 

Vexp -20% +48% +8% +40% 

Tkin +20% +13% +2% +22% 

+20% +13% +2% +22% 
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Fig. 6. Line profiles of the 845 - 752 transition for the various 
models indicated in Tables |7]and|8]for R Dor (top) and TX Cam 
(bottom). The left panels show the models with decreased pa- 
rameters; the right panels show models with increased param- 
eters. One parameter is varied, while the others are held fixed, 
thus leading to a change in the density profile when decreasing 
the expansion velocity. 



detailed abundance distribution is likely to have the largest ef- 
fect on the shape of the line profile. Detailed examination of the 
predicted HIFI line profiles based on the spectrally resolved 557 
GHz lines of R Cas, R Dor and W Hya, obtained with the Odin 
satellite, will be done in a forthcoming paper (Maercker et al., 
in prep.). Since this line is excited throughout the CSE, it will 
in particular help put better constraints on the size of the H2O 
emitting region. 

5.6. ISO co-lines 

The ISO LWS spectrum also contains high-lying transitions of 
the CO molecule (Fig.[TJ. These CO lines probe the same region 
as the H2O lines, while the radio lines probe the regions fur- 
ther out. Modelling these high excitation CO lines can be used 
as a consistency check for the adopted circumstellar model used 
in the ALI models. Due to weak lines and blending, in general 
only between one and three lines were found in the ISO spec- 
tra. W Hya, R Dor, and R Cas show the CO 7 = 16 ^ 15 and 
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Fig. 7. Eight ortho-HiO lines that can be observed with Herschel/HIFI TX Cam (top two rows), W Hya (middle two rows) and for 
IK Tau (bottom two rows) models based on the theoretical radius from Eq. |2] The dot-dashed line marks the stellar vlsr, the solid 
line shows vlsr ± Vexp- A number of line profile shapes are predicted, and the HIFI observations will help to verify or adjust the 
models. 
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y = 22— >21 lines. R Cas also shows the CO /= 17— > 16 line. 
In the remaining objects only the CO 7 = 16 — > 15 line could be 
found. For all objects except W Hya, the best-fit model fits the 
observations within +20%. The modelled CO lines for W Hya 
are ~ 80% weaker than the observations, in line with the higher 
mass-loss rate used by Justtanont et al. (12005b (see Sect. 13.2b . 

6. Discussion 

6.1. Comparison with previous results 

Abundance estimates for two of the sources treated here have 
been made previously. For W Hya, Justtanont et al. (12005b 
find an ortho-HaO abundance of 1.0 + 0.4 x 10"^ relative to H2. 
In their paper they use an ALI code to determine the H2O 
abundance by modelling of ISO SWS and LWS lines, and the 
Odin 557 GHz line. This agrees well with the abundance de- 
termined here (1.5 + 0.75 x 10"^). Our higher abundance is due 
to the lower mass-loss rate used (see Sect. 13.2b . Barlow et al. 
(11996b employ a large velocity gradient (LVG) model and de- 
termine a total H2O/H2 of 8 X 10 in the inner part of the CSE 
(r < 6 X lO''* cm) and 3x10""* in the outer part (r > 6 x 10'"* cm). 
Although the ortho-to-para ratio in their article is not very well 
constrained, between 1 and 3, their results are consistent with the 
abundances derived here. Zubko & Elitzur (12000b determine an 
abundance for ortho-H20 of 5.7 x 10"^, one order of magnitude 
lower than the results presented here, and a mass-loss rate one 
order of magnitude higher (2.3 x 10"^ M© yr"'). Their results are 
therefore consistent with ours, given the higher mass-loss rate. 
For R Cas, Truong-Bach et al. (119991 ) determine an ortho-H20 
abundance of only 7.9 x 10"^, two orders of magnitude lower 
than our result, using ISO LWS spectra. The mass-loss-rate used 
in their paper is somewhat lower than the one determined from 
the CO modelling here, further increasing the discrepancy in the 
abundance estimates. They do not model the dust emission ex- 
plicitly, but instead use a power law to reproduce the far-infrared 
continuum spectrum. They also estimate cooling by CO emis- 
sion by adopting a three level model, compared to the more 
complex scheme included here. On the other hand, they include 
a dust grain drift velocity (relative to the gas) that depends on 
the distance to the star Including a velocity field in our models 
would probably result in somewhat lower abundances. Although 
the models are different, this is far from enough to explain the 
different results. 

6.2. Evaporating pianetary systems? 

Assuming an ortho-to-para ratio of 1 gives an estimate of the 
total H2O abundance of ~ 6 x 10 in the majority of our stars. 
[Justtanont et al. (12005b and Zubko & Elitzur (12000b found an 
ortho-to-para H2O ratio of 1 and 0.8 for W Hya, respectively. 
Truong-Bach et al. ( 11999b derive a ratio of 3 for R Cas, in ac- 
cordance with what is expected in LTE.] The estimated W Hya 
H2O abundance is significantly higher than this, but, as outlined 
in Justtanont et al. (12005b . the W Hya mass-loss rate could be 
significantly higher than the value we use, and consequently the 
H2O abundance could be significantly lower The cosmic abun- 
dance of carbon and oxygen (Anders & Grevesse 11989b limits 
the abundance of H2O relative to H2 to not more than 1 x 10"^, 
assuming that all the carbon goes into CO, and all the remaining 
O goes into H2O. The abundance of H2O in thermal equilib- 
rium (TE) is approximately l-3x lO""* relative to H2 (Mamon et 
al. [T9B71 N ejad & Miller [T9881 Willacy & Millar [T9971 Duari 
et al. 11999b . Taken at face value, our results indicate that the 



abundances of H2O observed in M-type AGB stars cannot be ex- 
plained by TE chemistry. The uncertainty in our results is, how- 
ever, at least a factor of a few, and therefore does not rule out a 
TE abundance. NLTE chemistry can affect the H2O abundance, 
e.g., due to different densities and temperatures in shocked re- 
gions close to the star (Duari et al. 1 19991 Cherchneff l2006b . and 
photodissociation further out in the envelope will strongly af- 
fect the abundance (Mamon et al. 119871 Nejad & Miller , 19881 
Willacy & Millar [T997b . Other processes suggested to explain 
the relatively high water abundances are the evaporation of icy 
bodies (Justtanont et al. 120051) or grain surface reactions, such 
as Fischer- Tropsch catalysis on the surfaces of small metallic 
grains (Willacy i2004d . The former process was first suggested to 
explain the surprising detection of H2O in C-rich CSEs (Melnick 
et al. l200n Hasegawa et al. l2006b . 

Future observations with Herschel will make it possible to 
discriminate between the different formation scenarios, as sug- 
gested by e.g. Gonzalez-Alfonso et al. (12007b in the case of C- 
rich CSEs. A high abundance of H2O in the inner regions of the 
CSE will be observable through a significant amount of mid- 
excitation H2O transitions. These are already observed in the 
ISO spectra for our objects, indicating that at least some of the 
observed H2O is formed in the inner envelope. Both Fischer- 
Tropsch catalysis and the evaporation of cometary bodies result 
in similar shell-like H2O distributions. The evaporation hypoth- 
esis, however, implies the presence of HDO, the detection of 
which would rule out the formation of H2O from material ejected 
by the star 

6.3. The very iow abundance in WX Psc 

WX Psc has an estimated abundance that is two orders of mag- 
nitude lower than for the other objects. A possible explanation 
for this is depletion of H2O on dust grains. For a sample of car- 
bon and M-type stars, Schoier & Olofsson (2006) show that the 
abundance of SiO decreases in the CSEs with increasing Mj Vexp 
(a measure of the density), and explain this by depletion of SiO 
onto the dust grains. The freezing of water ice onto dust grains is 
also the explanation for the lack of H2O lines in the spectra of a 
sample of OH/IR stars (Justtanont et al. 12006b . Another possible 
explanation for the low abundance in WX Psc could be that the 
ISO water vapour observations are probing a present day mass- 
loss rate, whereas the CO transitions probe an earlier mass-loss 
epoch (Ramstedt et al. 12007b . If the present day mass loss is 
much lower, this would require an increase in the H2O abun- 
dance, thereby putting WX Psc more in line with the other ob- 
jects. Mass-loss-rate modulations in WX Psc with a low present 
day mass-loss-rate have been suggested previously by Kemper 
et al. (120031) and are further discussed by Ramstedt et al. (12007b . 

7. Conclusions 

We have demonstrated the ability of our accelerated lambda it- 
eration (ALI) code to model circumstellar H2O emission lines 
in M-type AGB stars at optical depths > 1000, and ortho-H20 
abundances are estimated. The models include the ground and 
first excited vibrational states, collisions within the ground vi- 
brational state, the first excited vibrational state, and between 
these states, and excitation by the stellar and dust radiation fields. 
The sample of stars modelled spans one order of magnitude in 
mass-loss rate. The ortho-H20 abundances determined, in the 
range (0.02-15)xl0"^ are estimated to be accurate to within 
50%, within the adopted circumstellar model, while the abso- 
lute error (including model-dependent assumptions) is within a 
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factor of 5. The estimated H2O abundances are high in compari- 
son to what can be produced by chemical models, but we explain 
that the estimates are rather uncertain. 

The outer radius of the H2O envelope is better determined 
for the low-mass-loss-rate objects and corresponds to within Icr 
to the radius given by the theoretical estimate in Eq. |2] 

Including excitation of the ground state through the vibra- 
tionally excited state affects the low-mass-loss-rate objects, in- 
creasing the model line fluxes by x 40%. This results in a lower 
abundance for these objects compared to models that do not in- 
clude the excited vibrational state (Justtanont et al. I2005I I. For 
higher mass-loss-rate objects the vibrational excitation is less 
important. 

Several formation processes for the observed H2O abun- 
dances are possible. LTE or NLTE models may explain the for- 
mation of H2O in the innermost layers of the CSE, whereas 
Fischer- Tropsch catalysis or the evaporation of Kuiper-belt like 
objects are possible sources of H2O in the outer parts of the enve- 
lope. Future observations with the Herschel telescope will help 
to diff'erentiate between these scenarios. The HlFl lines are gen- 
erally less saturated and will therefor give a better estimate of the 
circumstellar H2O abundance. In particular, optically thin lines, 
such as the 845 - 752 transition, are sensitive to the H2O abun- 
dance, and will help constrain this parameter. The ground state 
line at 557 GHz, on the other hand, is excited throughout the en- 
velope and will help in setting constraints on the size of the H2O 
envelope. 

For WX Psc, we suggest that sticking of H2O onto 
dust grains may explain the observed low amount of H2O. 
Alternatively, the present day mass-loss rate is low. 
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Appendix A: Accelerated lambda iteration (ALI) 

A.1. General considerations 

The radiative transfer models used in this work are created us- 
ing the accelerated lambda iteration (ALI) method (Rybicki & 
Hummer [19911 fT992b . As for most NLTE codes, the ALI tech- 
nique makes an initial guess of the level populations «, and 
solves the statistical equilibrium (SB) equations (for transition 



'YjniAn, - {ni,Bn - niBii,)Ju>] - 
i'<i 

^[«/'A/'/ - {niBu, - ni,Bri)Jw] + 
'^(riiCir -n,.Cn) = 0, 



/'>/ 



(A.l) 



where A//- is the Einstein A-coefficient, Bw and B^i are the 
Einstein B-coefficients for stimulated emission and absorption 
respectively, and Cw are the collisional rates. //// is the mean 
integrated intensity 



///' 



(A.2) 



where /(yu, v) is the specific intensity at frequencey v along direc- 
tion //, and (pii'ifi, v) is the normalized line shape profile (Doppler 
profile) used as a weight function. In the spherically symmetric 
case it is enough to solve the SE equations by determining 
at every radial point in the circumstellar shell, considering both 
the local and non-local contributions from all directions fi. This 
is done by introducing the lambda operator A^y, which relates to 
the specific intensity /(//, v) as 



I{fl, v) = A^v5tot(y", V) -t- Lbg, 



(A.3) 



where SiotiM' ^) is the total source function (including dust and 
all overlapping lines) and Li,g is the contribution from the back- 
ground to /(yU, v). The lambda operator can be seen as a N x N 
matrix with radial points. The diagonal elements of the matrix 
describe the local contribution to the specific intensity, while the 
off-diagonal elements describe the contribution from all other 
parts of the cloud. If t(jU, v) is the local optical depth in direction 
/i, the diagonal elements become 



Kr,r{H, V) = 1 - e 



(A.4) 



For high optical depths (t(jU, v) » 1) A^^ 1, whereas in the 
optical thin case {t{^, v) 0) A,- ,- — > 0. Instead of applying the 
lambda iteration directly, the accelerated lambda iteration tech- 
nique uses an approximate lambda operator A*y to take care of 
the high optical depths separately 



/(yU, v) = A*y5 , 



,(yU, v) + (A^„ - A* )5;„,(yu, v). 



(A.5) 



5'tot ' (A', v) is the source function from the previous iteration and 
can be used, since convergence guarantees that 5 tot ' — * Stot- It is 
convenient to simply choose the diagonal elements of the lambda 
operator as the approximate lambda operator This makes the 
matrix inversion easily possible, and the SE equations can be 
solved at every radial point individually, based on the approx- 
imate lambda operator only. Using the accelerated lambda it- 
eration technique improves the rate of convergence compared 
to normal lambda iteration, particularly in cases of high op- 
tical depths. To further improve convergence, Ng acceleration 
(Ng |1974] i is used, with a good description of the method given 
by Olson et al. (fW86l l. 
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Fig. A.l. Results from benchmarking our ALI code against the 
results found in van Zadelhoff et al. ( I2002l i for their problems 2a 
and 2b, consisting of an inside-out collapsing spherical cloud. 
The figure shows the radial distributions of the relative num- 
ber density of HCO^ molecules in the rotational levels 7=1 
and 7 = 4 as obtained by our ALI code after convergence. Also 
shown are the deviations (in percent) of our results when com- 
pared to the mean of the results found by all other radiative trans- 
fer codes, as presented in van Zadelhoff et al. ( I2002l i. 



A.2. Benclimarking 

When developing a new radiative transfer code it is impor- 
tant to perform detailed tests as to its accuracy. A num- 
ber of such test problems are available for download at 
www.strw.leidenuniv.nl/~radtrans. Our ALI code has 
been benchmarked against the solutions to these problems as 
found by seven different radiative transfer codes (presented in 
van Zadelhoff et al. 12002b . We find that our ALI code repro- 
duces the results for all these test problems with an accuracy 
that matches all the other codes included in van Zadelhoff et al. 
(120021) . 

As an example. Fig. lA. II shows the solution obtained from 
our ALI code for problems 2a and 2b of van Zadelhoff et 
al. ( I2002I I. These problems consist of an inside-out collapsing 
spherical cloud and contain radial gradients in most physical 
parameters such as velocity, density, temperature and micro- 
turbulent velocity. The excitation analysis is performed for 
HCO^, and the final level populations (in relation to the total 
number density of HCO^ molecules, «tot) for the 7=1 and 7 = 4 
rotational levels are shown in Fig. I A. II Also shown are the de- 
viations of our ALI code (in percent) from the mean of the re- 
sults obtained by the other radiative transfer codes presented in 
van Zadelhoff et al. ( 120021 ). The spread in the results, around the 
mean, between various codes in van Zadelhoff et al. ( 120021 ) is for 
the 7=1 level about ± 1 % and + 2% for problems 2a and 2b, re- 
spectively. For the 7 = 4 level a larger spread of about +5% and 
+ 15% was obtained for problems 2a and 2b, respectively. 



